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Abstract 

We demonstrate an application of spherical harmonic decomposition to anal- 
ysis of the human electroencephalogram (EEG) . We implement two methods 
and discuss issues specific to analysis of hemispherical, irregularly sampled 
data. Performance of the methods and spatial sampling requirements are 
quantified using simulated data. The analysis is applied to experimental EEG 
data, confirming earlier reports of an approximate frequency-wavenumber re- 
lationship in some bands. 
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Typeset using REVTeK 



I. INTRODUCTION 



The human electroencephalogram (EEG), as measured at the scalp, represents a super- 
position of electric fields resulting from post-synaptic potentials in neocortex, the thin (2 to 
5 mm) surface layer of human brains. Several models of neocortical dynamics treat EEG as 
a mixed global/local phenomenon and a better understanding of its spatial-temporal 

dynamics is necessary for evaluation and refinement of these models. Its temporal behavior 
has been studied at length, both by clinical observation and with such tools as power 
spectra coherence [^, the Hilbert transform [0, and many others. However, until re- 
cently poor spatial resolution (due to minimal electrode sampling and under-use of head 
models) has limited spatial analysis of EEG 

As a potential field on a near-hemispherical surface, EEG is amenable to analysis by 
spherical harmonic decomposition. In this paper, we apply two methods of decomposition 
(one described by Cadusch and one adapted from Misner |Ty]) to 131-channel EEG data. 
Using simulated data, we discuss issues and pitfalls relevant to such an analysis, specifically 
the effects of limited and irregular sampling density, integration over a hemisphere, and 
deviations from a spherical surface. From the experimental data, we then draw conclusions 
regarding the frequency-wavenumber relation of neocortical activity. 

II. METHOD 



We use the real spherical harmonics [jTT|, defined on the sphere f2 and described by the 
orthogonality integral 



In theory, a potential field $(1^) may be decomposed into spherical harmonic amplitudes 
^im defined by 

^im = J YU0AMdA)d^^- (2) 
n 

In the example of EEG and similar data we encounter three major and two minor issues. 



A. Sampling 



First, when attempting to decompose experimental data, we sample $(f^) at specific 
locations F. Assuming near-regularly spaced electrodes, our maximum resolvable / is deter- 
mined by a spherical analog of the familiar Nyquist limit [|12| fmax = 1/(2 AT). With mean 
angular inter-electrode distance 7, we initially adopt a conservative limit of 
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(3) 



or Imax = 6 for our 131-channel electrode cap. Analog pre-filtering to avoid spatial aliasing 
is not required here due to the low-pass characteristics of the head volume conductor |13|. 
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B. Regularization 



With sampled data F{r), the discretized version of the decomposition in Eq. ^is unstable 
in higher Z-indices. (An apparently accurate reconstruction of the signal may be generated, 
with large artifacts in the higher spatial frequencies.) We must invoke constraint or regular- 
ization techniques to address this issue. Cadusch et al approached the problem as a side 
issue of spherical spline interpolation. The estimate is constrained by the spline con- 
straints, and the problem for a given sampling grid and Imax is reduced to a multiplication 
by a matrix of fiim coefficients: 

^'im = ^/i«m(x)F(x). (4) 



Recently, Misner |T^ introduced a more complex method for decomposition on a rect- 



angular three-dimensional grid; that is, generalized to use Ynim{r,9,(f)). We implement here 
the special case of sampling on a spherical surface, more relevant to EEG analysis, in which 
r is constant. In Misner's method, a correction matrix of Gab accounts for discretization 
and limited sampling: 

Gab = Yu yA{x)YBix)w^, (5) 

where A and B refer to index groups (Im). Here we use the real harmonics and replace 
Misner's weight function with the effective area of each electrode. G"^^, the matrix 
inverse of the Gab, is used to calculate "adjoint spherical harmonics" Y^. Finally, a set of 
coefficients Rim, analogous to Cadusch' /i;^, is generated and used as in Eq. | to estimate 



C. Hemispherical sampling 

Particularly relevant to EEG analysis is the error introduced by sampling over only half 
of the sphere. This sampling corresponds to only ^ of a spatial cycle of the / = 1, m = 
function, suggesting potential accuracy problems for functions involving / = and 1. It is 
also clear that the functions Yim will no longer be orthonormal for < 6 < 7t/2 only; rather, 
we replace the 6 in Eq. ^ with an error quantity e: 



2lT 2 








In general, our hemispherical estimates will be related to the hypothetical full-sphere 
result by the matrix of e^^'s: 

$' = e<|) (7) 

and our become somewhat ambiguous between certain sets of (/,m). Although it 
may seem appropriate to invert e and calculate a more accurate result, the matrix is ill- 
conditioned {R > 10^, where the 2-norm condition number R is the ratio of the largest 
singular value of e to the smallest) and thus the inversion is problematic. 
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In addition, if we use the hemispherical region to calculate Misner's Gab in Eq- @, the 
resulting matrix of Gab is ill-conditioned [R > 10^) and thus the G^^ cannot be reliably 
found. Rather, we created a mirrored set of electrodes F', calculated the matrix R of Rim 
for sample set (T U V), and discarded the antipodal rows of R. Cadusch' spline method, 
while still subject to Eq. |^ is native to the hemispherical surface and requires no further 
manipulation. 



D. Coordinate orientation 



In many problems, the sphere has no preferred direction. The m-indices are usually col- 
lapsed to produce an angular power spectrum estimate G (/) as a function of wavenumber 



/: 

G{1) = Y: {^'imf (8) 
m=—l 

which is independent of coordinate system orientation. As well, we found the "hemispherical 
error" in /-spectrum to be independent of orientation. In some EEG studies, of course, the 
orientation of the underlying cerebral hemispheres may be relevant. In such cases, local 
spatial Fourier analysis |1| should adequately complement our decomposition without the 
complication of distinguishing m-modes. 



E. Non-spherical media 

We assume that our medium is a sphere, whereas biological data is often sampled on an 



irregular surface. The upper surface of the "average" human head [T^ may be represented 
as a hemiellipsoid with axes a = 10.52 cm, b = 7.66 cm, and c = 8.41 cm, or alternatively 
25%, -9%, and 0% elongation from a perfect sphere. Although prolate spheroidal harmonics 



have been applied to biophysical field problems [|T6|-[T8[], the technique is often unwieldy. In 



comparison to error from e^s, especially for low /, we assume the error due to approximating 
the ellipsoidal surface with a spherical surface is negligible. 



III. APPLICATION TO SIMULATED DATA 

We generated evenly tessellated, hemispherical electrode maps of 74, 187, 282, and 559 
electrodes, in addition to common experimental maps of 20, 64, and 131 electrodes. Five 
hundred potential maps were simulated for each electrode configuration. Each potential 
map was randomly generated with harmonics of degree / = 6, such that the varied with 
uniform distribution between and 1. Power spectrum estimates G (/) were then calculated 
for each map, using both methods. Figure |T] shows Pearson's correlation coefficients ri, 
calculated between G (/) and G (/) over the 500 trials for each electrode map. 

We have noted that the error due to cab causes power from one (/, m) -component to be 
misinterpreted as power in another, often of different /. Therefore, we might expect either 
method's performance to depend on the /-spectrum being analyzed. Using preliminary 
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experimental data, we constructed an approximate power spectrum Gnormi}) for average- 
referenced scalp EEG, peaked at / = 1 and / = 2, and decaying with thereafter. Another 
five hundred potential maps were generated, with uniformly distributed between 

< < %^ (9) 

to simulate a physiologically realistic distribution of /-spectra. Power spectrum estimates 
G (/) were calculated for each map using both methods. Results are shown in Fig. |^. 

In general, results for the spline method — though often quite accurate — were dependent 
on the distribution of /-spectra being measured, exact electrode positions, and electrode 
numbers. Results for the adjoint harmonic method seemed more robust, even for sparse 
(n = 64) sampling, although accuracy was somewhat less in the higher harmonics. 

In both methods, for / = 6 we observed minimal improvement for more than 131 elec- 
trodes. We thus believe that our 131-channel sampling is an appropriate tool for further 
study. Furthermore, given the limit in Eq. ^ and the known volume-conductor attenuation 



of higher modes |]I3| , we suggest that study of spatial frequencies higher than approximately 
/ = 8 will be better served by intracranial EEG than by denser electrode maps. 

In general, for low / the adjoint harmonic method seemed more consistent. We examined 
typical 131-channel decompositions (Fig. ^) to investigate further. Both methods accurately 
reproduced the potential maps (r > 0.9 for 131 channels). The spline method, however, 
was slightly unstable for low /, and the erroneous negative are reflected in the power 
spectrum. 



IV. REFINEMENTS AND ANALYSIS OF ERROR 



Any application of the spherical harmonic decomposition should take into account the 
estimated relative contribution of various error sources. Aside from measurement and ex- 
perimental error, these may be divided into three categories: sampling error, orientation 
error, and hemispherical error. 



A. Sampling error 



Figures |I| and |^ indicate minimal improvement for / < 6 with more than 131 electrodes. 
We can thus deduce that the Nyquist-like limit in Eq. ^ is an appropriate guideline. When 
using coarser sampling, we expect some decrease in performance for higher /. Decreased 
accuracy for 20, 64, and 74 channels (seen in Figs. |l] and ^, particularly for 20 channels) 
may be attributed to sampling error. 



B. Orientation error 



For a given /-spectrum, results will vary if power is randomly distributed across the m's; 
that is, various m-components interact differently with our hemispherical sample grid. Five 
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hundred random /-spectra, with reahstic distribution of G{1), were generated. Thirty 131- 
channel maps, randomly varying in m-power, were generated for each original /-spectrum. 
Figure ^ displays the resulting accuracy, again shown as correlation coefficient ri between 
G{1) and G (/) over the 100 trials. 

Very little is gained in this simulation by decomposition of multiple epochs. As de- 
scribed in the following section, error due to m-distribution of power is largely swamped 
by hemispherical error. In practice, however, we must emphasize (in the presence of ran- 
dom measurement noise) the importance of averaging decompositions across many epochs. 
Orientation error will also become significant if our sampling grid is severely non-uniform. 

C. Hemispherical error 

In Sec. H above, we have discussed the hemispherical error cab- Although it is impos- 
sible to improve our decomposition results by inverting the matrix e, we may generate a 
corresponding matrix for the power spectrum result and use it to estimate the contribution 
of hemispherical error. 

Power in a single harmonic ¥[^{0, <p) is blurred by the hemispherical decomposition into 
surrounding harmonics. Using the 131-channel sampling map, we generated five hundred 
potential maps for each of / = ... 6, each with one unit power distributed randomly among 
the available m's. By averaging over the five hundred resulting power spectra, for each /, 
we obtained an empirical "averaged blurring matrix" E for power spectra obtained from 
hemispherical decomposition. That is, 

G(/)^EG(/). (10) 

The typical E for both methods is a blurred identity matrix; that is, error in power spectra 
is largely between adjacent /. It is again tempting to invert E, de-blur our spectra, and cal- 
culate a more accurate result, but although most E are invertable we found that for realistic 
spectra the benefit was marginal at best. Instead, E may be used to better understand the 
implications of hemispherical error. 

We calculated correlation coefficients as in Fig. ^, between EG(/) and G (/) over the 
500 trials of 30 epochs for 131 electrodes. The resulting higher correlations (although not 
applicable to a decomposition of real data) are plotted in Fig. ^ By comparison with Figs. 
^ and ^, the result indicates the importance of hemispherical error. In particular, after 
examination of typical E and e matrices, we may interpret the decreased performance at 
low / as blurring between adjacent wavenumbers. Furthermore, the increased effect, seen 
in Fig. 1^, of averaging across various m-distributions indicates that some abrupt changes 
in performance may be attributed to sensitive interactions between E-blurring and random 
m-distribution. 

Practically, the near-identity character of E is extremely useful. Hemispherical error 
manifests as blurring between adjacent /. Thus, we may expect composite measures such 
as the sum of power in /=0,1 to be substantially more accurate than individual estimates. 
Figure ^ displays the accuracy of /=0,1 and /=2. . . 6 adjoint harmonic power estimates 
(used below in our experimental trials), for 500 epochs, realistic /-distribution, and various 
sampling densities. 
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V. APPLICATION TO EXPERIMENTAL DATA 



Nunez in 1974 [0,^] and Shaw in 1991 [^], using Fourier analysis along linear electrode 
EEG arrays, observed a relationship between increasing spatial and increasing temporal 
frequency in the 8-13 Hz band, roughly consistent with simple wave dispersion relations. 
We attempted to duplicate this result in order to test the adjoint harmonic method under 
experimental conditions. We analyzed 131-channel EEG (resting, eyes closed) in five human 
subjects. Temporal Fourier coefficients were determined for 300 to 600 one-second epochs 
(depending on available data), and /-spectra averaged over these epochs. 

Results are summarized in Fig. |^ as the ratio of power in low (/ = 0, 1) to power in high 
(/ = 2,3,4,5,6) spatial frequencies. Above approximately / = 8 Hz, with increasing / we 
observed a general trend towards power in higher /. We also observed high-wavenumber 
power in the delta band (/ < 3 Hz). The alpha band (c. 8-13 Hz) was characterized by the 
highest power in low spatial frequencies. 

In order to rule out methodological artifact, we generated and analyzed 300 seconds of 
simulated EEG using 3602 uncorrelated sources, each generating 1/f noise through a head 
volume conductor model as described in [1^. As expected, the EEG-like noise (labeled RND 
in the figure) showed no relation between spatial and temporal frequencies. 



VI. CONFIDENCE INTERVAL ESTIMATION 

Estimates of the temporal power spectrum are known to vary in chi-square distribution 
T^, assuming normally distributed estimates of the underlying Fourier coefficients. Error 



distribution for spatial spectrum estimates, on the other hand, is complicated by the de- 
pendence of hemispherical and orientation error on the entire /-spectrum. For composite 
measures of both spatial and temporal spectra, such as shown in Fig. |^, the situation be- 
comes even more problematic. We propose an empirical test for estimation of such confidence 
levels, analogous to the randomization tests commonly applied in nonparametric statistical 



analysis |]22 
Let 

A = P^, (11) 

where Gqi is the total power in harmonics /=0 and 1=1, and 6*2. ..6 is the total power in 
harmonics 1=2 through 1=6. Let Gqi, G2...6y and A represent estimates of the same. Above, 
we calculated Aj for various temporal frequencies. Here, we calculate an approximate 95% 
confidence interval for single-epoch estimates of the actual Af. The confidence interval will 
apply only to the spatial spectrum composite measure, neglecting error (or nonstationarity) 
in temporal frequency spectra, which for many applications may be as important. Note, 
though, that for 300 epochs the normalized standard error of a temporal power spectrum 
estimate is less than 6%. 

To determine an empirical confidence interval, we would typically examine the distri- 
bution of random re-samples. In this application, we created many random /-spectra from 
an estimated distribution of /-power, simulated many decompositions, and examined the 
resulting distribution as follows. 
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Since hemispherical error is dependent on /-spectrum, the result will be influenced by 



the distribution of the random spectra. Srinivasan et al [|13| analytically estimated the 
spatial frequency domain transfer function for volume-conduction blurring of scalp potential 
as proportional to (2/ + l)~^. This "spatial smearing" is due mainly to the poorly conducting 
skull and physical separation between cortical current sources and scalp electrodes. In our 
calculation, we assumed that underlying /-spectra vary in uniform distribution in proportion 
to (2/ + and that with average- referenced data the contribution of /=0 is negligible. A 
large number (20,000) of /-spectra were generated, randomly selecting for each /-bin a value 
from the appropriate distribution. The decomposition was performed, and the composite 
measure A calculated, for each randomized spectrum. 

By examining the distribution of known surrogate Arand which produce a certain estimate 
Arand , wc cau estimate an empirical confidence interval for our spectral estimate. In Fig. |^, 
we show the scatter plot of Arand against Arand with 95% confidence intervals. For a given 
estimate A and the assumptions discussed above, 95% of the time, the actual A will fall 
between the two lines shown. 

A similar procedure may be used to calculate confidence intervals for other measures, 
whether the actual Gi or other composite measures. Careful judgment must be applied 
when estimating confidence intervals for multiple-epoch measures such as shown in Fig. |^. 
As demonstrated earlier in this paper, variation in the m-component of an /,m-spectrum 
only allows us to "average out" the minimal orientation error. Variation in hemispherical 
error (dependent on /-spectrum), without gross violation of the stationarity assumption, is 
necessary for the average of estimates A over multiple epochs to converge to A. 

VII. DISCUSSION 

Our simulations provide a firm basis for application of spherical harmonic decomposition 
to irregularly sampled, hemispherical data such as EEG. Our hemispherical modification 
of Misner's adjoint harmonic method [1^ proved most consistent. However, for physiologi- 



cal data of known power distribution, the spline method is complementary and may be 
slightly more accurate with high-density sampling. It seems that, within the conservative 



band-limit of equation 0] and the known spatial filter properties of the head |]T3[, decom- 
position accuracy will not be materially improved by more than 131 electrodes for scalp 
EEG. We suggest that confidence intervals for such decompositions, or for decomposition- 
derived measures, be determined empirically using randomized data. Furthermore, while 
single-decomposition errors are relatively large, with multiple epochs the experimental ac- 
curacy may be increased substantially. For this averaging to be both valid and effective, we 
must assume a quasi-stationary wavenumber spectrum across our epochs, but with sufficient 
random variation in hemispherical error for our estimates to converge upon the mean. In 
addition, especially in EEG applications, we must remain aware of the limitations inherent 
in collapse across m's (we assume the orientation of the underlying cerebral hemispheres is 
irrelevant) and the use of spherical harmonics on a hemispheroidal surface. 

The dynamical properties of human EEG rhythms are quite complicated, varying sub- 
stantially between individuals and brain states. Furthermore, physiologically-based theo- 
retical models point to substantial nonlinear effects and interactions across spatial scales 
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Pj23|-p6|. Despite all the obvious complications, results from the spherical harmonic de- 



composition of experimental EEG agreed qualitatively with crude linear electrode array 
results |1|J2^. These results were seen in all subjects and are consistent with a mixed 
global/local model of cortical dynamics, in which lower global mode oscillations produce 
alpha rhythm, superimposed on local (spatially uncorrelated) activity in various frequency 
bands 0. Further study of spatiotemporal EEG dynamics, using spherical harmonic de- 
composition, should shed more light on these issues. 
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FIGURES 



FIG. 1. Correlations between actual and estimated /-power for uniformly distributed random 
spectra. Five hundred potential maps were generated from known wavenumber spectra, with 
random power in each /-component, uniformly distributed between and 1. For each of seven 
electrode densities, wavenumber spectra were estimated by spherical harmonic decomposition of 
the 500 sampled maps. Correlations between actual power and estimated power were calculated 
over the 500 trials for each /-component. Shown here for (a) adjoint harmonic and (b) spline 
methods, these correlations are a measure of the quality of a single decomposed power spectrum. 

FIG. 2. Correlations between actual and estimated /-power for more realistically distributed 
random spectra based on genuine EEG data. As in Fig. |, but in original spectra random power 
in each /-component is uniformly distributed between and (2/ + 

FIG. 3. Topography (left column), /,m-spectra (center column), and /-power (right column) 
for a typical 131-channel spherical harmonic decomposition. The original map is shown in (a). The 
adjoint harmonic method (b) reconstructs topography and gives an approximation of /-spectrum. 
Although the spline method (c) also reconstructs potential topography, we observe irregularities in 
the lower / amplitude estimates that contribute to decreased performance for these wavenumbers, 
and a less accurate /-spectrum estimate. 

FIG. 4. Correlations between actual and estimated /-power for multiple-epoch, ad- 
joint-harmonic estimates of the same /-spectrum, with epochs varying only in m-component. For 
a reasonably isotropic and dense sample array, such as the 131-channel EEG grid used here, there 
is little orientation error and thus little improvement in results. 

FIG. 5. Sampling a full sphere with 262 channels and the adjoint harmonic method, correlations 
between actual and estimated /-power are shown for multiple-epoch (varying only in m-component) 
estimates of the same /-spectrum. By sampling over the full sphere, we eliminate hemispherical 
errors illustrated in Fig. ^ Remaining errors are due to orientation (note improvement with 
multiple epochs) and imperfect sampling. 

FIG. 6. Correlation coefficients, using the adjoint harmonic method, obtained by comparisons 
of estimated to actual summed power measures. The solid line represents power in /=0 and 1=1 
modes, and the broken line represents power summed over modes 1=2 through 1=6. Increased 
accuracy (as compared to part A of Fig. ^) is because most hemispherical error manifests as 
blurring between power in adjacent Vs. 
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FIG. 7. Five to ten minutes of resting, eyes-closed EEG were collected with 131 channels from 
each of six subjects (one duplicated). Complex temporal Fourier coefficients were calculated for 
one-second epochs and subjected to spherical harmonic spatial decomposition using the adjoint 
harmonic method. Resulting wavenumber spectra were averaged for each 1-Hz band over the 300 
to 600 epochs. The ratio of power in Z=0,1 to Z=2,3,4,5,6 is plotted as a simple indicator of a bias 
toward higher spatial frequencies at higher temporal frequencies (greater than about 10 Hz). This 
result is qualitatively consistent with the postulated existence of an approximate EEG dispersion 
relation, perhaps with alpha rhythm (8-13 Hz) representing the fundamental and lower overtones. 
A surrogate signal (dotted line), composed of random EEG-like noise and subjected to the same 
analysis, showed no such relation. 

FIG. 8. 95% confidence intervals for single-epoch estimates of the power ratio 
A = G'/=o,i/G'/=2...6- Twenty thousand potential maps were generated from known, random, real- 
istically distributed (based on genuine EEG data) /-spectra, and decomposed using 131 channels 
and the adjoint harmonic method. Here, known A are plotted against the resulting estimated 
A. Solid lines indicate the empirical 95% confidence interval for a given estimate of A. Multi- 
ple-epoch estimates will result in much smaller intervals, depending on the variation in Z-spectra 
being decomposed. 
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